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We formulate a quantum Monte Carlo (QMC) method for calculating the ground state of many-boson systems. 
The method is based on a field-theoretical approach, and is closely related to existing fermion auxiliary-field 
QMC methods which are applied in several fields of physics. The ground-state projection is implemented as a 
branching random walk in the space of permanents consisting of identical single-particle orbitals. Any single- 
particle basis can be used, and the method is in principle exact. We illustrate this method with a trapped atomic 
boson gas, where the atoms interact via an attractive or repulsive contact two-body potential. We choose as the 
single-particle basis a real-space grid. We compare with exact results in small systems, and arbitrarily-sized sys- 
tems of untrapped bosons with attractive interactions in one dimension, where analytical solutions exist. We also 
compare with the corresponding Gross-Pitaevskii (GP) mean-field calculations for trapped atoms, and discuss 
the close formal relation between our method and the GP approach. Our method provides a way to system- 
atically improve upon GP while using the same framework, capturing interaction and correlation effects with 
a stochastic, coherent ensemble of non-interacting solutions. We discuss various algorithmic issues, including 
importance sampling and the back-propagation technique for computing observables, and illustrate them with 
numerical studies. We show results for systems with up to N ~ 400 bosons. 



I. INTRODUCTION 

The study of many-body quantum systems has been a very 
challenging research field for many years. Computational 
methods have often been the way of choice to extract theoreti- 
cal understanding on such systems. Most computational quan- 
tum mechanical studies are based on simpler mean-field theo- 
ries such as the Gross-Pitaevskii (GP) equation for bosons or 
the Kohn-Sham density-functional theory (DFT) for fermions. 
Despite their remarkable success, the treatment of particle 
interaction or correlation effects is only approximate within 
these approaches, and can lead to incorrect results, especially 
as the strength of particle interactions is increased. It is there- 
fore necessary to develop alternative computational methods 
that can describe the effect of interaction more accurately and 
reliably. 

In this paper we present a quantum Monte Carlo (QMC) 
method to study the ground state of many-boson systems. 
The method is in principle exact. Our interest in the devel- 
opment and use of this method was motivated by the real- 
ization of the Bose-Einstein condensation in ultracold atomic 
gases 1 1]. These are dilute gases consisting of interacting al- 
kali atoms. The interaction among the atoms is well described 
by a simple two-body potential, either attractive or repulsive, 
based on the scattering length. For weakly-interacting sys- 
tems the mean-field GP approach has, as expected, performed 
extremely well OB]. More recently, Fesbach resonances 0] 
have successfully been used as a powerful way to tune the 
strength of the interaction experimentally. This provides a 
source of rich physics, and increases the need for theoretical 
methods which can benchmark GP and provide an alternative 
where GP is inadequate. 

Several QMC methods exist for calculating the properties 
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of interacting many-body systems. The ground-state diffu- 
sion Monte Carlo |5 1 and the finite-temperature path-integral 
Monte Carlo (PIMC) 0] methods, which work in many- 
particle configuration space and in the first-quantized frame- 
work, have been successfully applied to a variety of boson and 
fermion systems. In the context of atomic gases, Krauth 01, 
Gruter et al. H, and Holzmann and Krauth ||9|] have em- 
ployed PIMC to study finite-temperature properties of trapped 
bosons with positive scattering lengths, modeling the two- 
body interactions by a hard-sphere potential. Glyde and co- 
workers have studied the ground state of trapped bosons, also 
by hard spheres fTolfTlll . Ulmke and Scalletar fl2ll did finite- 
temperature QMC calculations on quantum spin systems and 
the Bose-Hubbard model. In the latter calculation, a hard-core 
repulsive potential was assumed, which allowed a transforma- 
tion of the problem into an XXZ spin-like problem that can be 
treated with a fermion QMC method. 

Our method is based on the auxiliary field quantum Monte 
Carlo (AFQMC) approach 0H. The AFQMC is a field- 
theoretical method, where many-body propagators resulting 
from two-body interactions are transformed, by use of auxil- 
iary fields, into a many-dimensional integral over one-body 
propagators The many-dimensional integral is 

then computed using stochastic means. The AFQMC frame- 
work is appealing for several reasons. Working in second- 
quantization, it automatically imposes the proper particle- 
permutation symmetry or antisymmetry. It provides a many- 
body method with close formal relation to mean-field ap- 
proaches, as we discuss later. In addition, it allows convenient 
calculation of the observables and correlation functions. 

The AFQMC method has been widely employed to study 
fermion systems in condensed matter 1171 [Ha U^l . nuclear 
physics 12(1 12111 . and lattice gauge theory. In this paper, we 
generalize the fermion ground-state auxiliary-field quantum 
Monte Carlo method fl9l l22ll to many-boson systems. We 
project the many-body boson ground-state from an initial trial 
state I ^> T ) . Our choice of | ^f T ) is a permanent consisting of A*" 
identical single-particle orbitals, which was first suggested in 
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a model calculation by Sugiyama and Koonin 11411 . The many- 
body ground state is projected from \^ T ) with open-ended, 
branching random walks to sample the auxiliary fields. We 
formulate an importance sampling scheme, which greatly im- 
proves the efficiency of the method and makes possible sim- 
ulations of large systems. We also discuss in detail the back- 
propagation technique which allows convenient calculation of 
virtually any ground-state observables. 

Our method retains all the advantages of AFQMC. It al- 
lows the use of any single-particle basis, which in this paper 
is chosen to be a real-space grid. As we discuss in Sec. lVIl it 
provides a means for true many-body calculations in a frame- 
work which closely relates to the GP approach. The approach 
can be viewed as a stochastic collection of parallel GP-like 
calculations whose "coherent" linear combination gives the 
interaction and correlation effects. 

In this paper we present our QMC method for bosons and 
discuss its behavior and characteristics. We use a trapped 
atomic boson gas as our test system, where the atoms inter- 
act via an attractive or repulsive contact two-body potential. 
A sufficiently detailed description of the method is given to 
facilitate implementation. Compared to its fermionic counter- 
part, our method here is formally simpler. It therefore also 
offers opportunities to study algorithmic issues. Because of 
the intense interest in methods for treating correlated systems 
(fermions or bosons) and the relatively early stage of this type 
of QMC methods, a second purpose of the paper is to use 
the bosonic test ground to explore, discuss, and illustrate the 
generic features of ground-state QMC methods based on aux- 
iliary fields. An example is the case of repulsive interactions, 
where a phase problem appears in a bosonic system, which 
provides a clean test ground to study methods for controlling 
this problem ll22ll . which is crucial for applications in fermion 
systems. The majority of the applications in this paper will be 
to systems where exact results are available for benchmark. 
These include small systems, which can be diagonalized ex- 
actly, and the case of untrapped bosons with attractive inter- 
actions in one dimension, where analytical solutions exist. It 
is worth emphasizing that the method scales gracefully (sim- 
ilar to GP) and allows calculations for a large number (N) 
of bosons. We will show results for larger systems (~ 1000 
sites and hundreds of particles) in one- and three-dimensions 
to illustrate this. 

Our paper is organized as follows. In section |n] we es- 
tablish some conventions and review the basic ground-state 
projection and auxiliary-field quantum Monte Carlo method. 
In section [TO] we introduce our new AFQMC implementa- 
tion for bosons, including the formulation of an importance- 
sampling scheme and the back-propagation technique for con- 
venient calculation of virtually any ground-state observables. 
In section Hvl we describe the implementation of our method 
to study the ground state of a trapped Bose atomic gas, which 
we model by by a Bose-Hubbard Hamiltonian with an exter- 
nal trapping potential. We also describe our implementation 
of the GP approach to study the same Hamiltonian. In section 
fvl we present our computational results. We benchmark the 
method in systems where exact results are available. We also 
provide examples to illustrate the behavior and key charac- 



teristics of our method. We carry out GP calculations on the 
same Hamiltonian and compare the results with those from 
our QMC calculations. In section fVll we comment on some 
characteristics of the method, further discuss its relation to 
and differences from GP, and mention future directions and 
some immediate applications of this method. Some comput- 
ing issues will also be discussed. Finally, in the appendices 
we provide additional technical details of the method. 



II. BACKGROUND 
A. Many-body Hamiltonian 

We use the second quantized formalism throughout this pa- 
per. We assume that an appropriate set of single-particle basis 
{IXi)} nas been chosen, in terms of which the wave func- 
tions will be expanded. For simplicity, we assume that the 
single-particle basis is orthonormal, although this is not re- 
quired. The number of basis states is M. The operators cj and 
Cj, respectively, are the usual creation and annihilation oper- 
ator for the state They satisfy the commutation relation 
[Cj,ct]_ = 5^. This automatically imposes the symmetriza- 
tion requirement of the many-body wave functions. 

We limit our discussion to a quantum-mechanical, many- 
body system with two-body interactions. The Hamiltonian H 
has a general form of 

H = K + V, (1) 

where K is the sum total of all the one-body operators (the 
kinetic energy and external potential energy), 

K A', , (•; ,•, : 

and V contains the two-body interactions: 

ijkl 

Our objective is to calculate the ground state properties of 
such a system, which contains a fixed number of particles, 
N. 



B. Ground state projection 

The ground state wave function |$ ) can be readily ex- 
tracted from a given trial solution \^> T ) using the ground-state 
projection operator 

V gs = e - Arfl e ATE ^ , (2) 

where Et is the best guess of the ground-state energy, pro- 
vided that \^ T ) is not orthogonal to |$ ). Applying the oper- 
ator Vgs repeatedly to the initial wave function would 
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exponentially attenuate the excited-state components of the 
initial wave function, leaving only the ground state: 



where v is a one-body operator: 



(Pgs) n |*T> 

Wo) 



l^o) 
l*o) 



(3a) 
(3b) 



Because of its resemblance to the real-time propagator, the 
operator "Pgs is also called the imaginary-time propagator. In 
ground-state QMC methods, T^gs is evaluated by means of a 
Monte Carlo sampling, resulting in a stochastic representation 
of the ground-state wave function. 



C. Basic auxiliary-field method 



The hermiticity of V allows us to decompose it into a sum 
of th e sq uare of one-body operators {fi,} (see, for example, 
Refs.EaandEl: 



V 



(6) 



Two essential ingredients are needed in order to evaluate 
Pgs within a reasonable computin g ti me. The first is the 



Trotter-Suzuki approximation 12 31 12411 . The propagator is 
broken up into a product of exponential operators, which be- 
comes exact in the limit At — ► 0. The second-order form of 
this approximation is 



-At(K+V) 



e -|ArK e -ArV e -iArir 



0{At 3 ). 



(4) 



The second ingredient is the Hubbard-Stratonovich (HS) 
transformation fl5l flfill . which allows us to reduce the two- 
body propagator to a multidimensional integral involving only 
one-body operators, using the following identity: l25ll 



;Atv 



1 



'2-k J- 



dx e 



x\/ At v 



(5) 



Because of this, we can always apply the Hubbard- 
Stratonovich transformation on a general two-body potential 
operator: 



n 



e 2 



_ 1 2 

e 



(7) 



TT / da^l— e *^«* + 0(At 2 ) . 



In general, the Trotter breakup incurs an additional systematic 
error of C(Ar 2 ). 

Applying these two procedures, we obtain an approximate 
expression of the ground-state projection operator: 



,AtB t 



,-iArJf 



n 



AriC 



0{At 2 ). 



(8) 



where p(x) is the normalized Gaussian probability density 
function with unit standard deviation: p(x) = —)=e~2 x . 
This approach is applicable to both boson and fermion sys- 
tems. It enables us to compute the exact ground state of a 
quantum many-body system. To reduce the systematic er- 
ror from the finite timestep At, the so-called "Trotter er- 
ror", small timesteps At are necessary. Often, calculations 
are performed for several At values, then an extrapolation to 
Ar — ^ is made to remove the Trotter error. 

For convenience we define the following notations: 

• x = {x l ,x 2 , . . .} : the collection of all the auxiliary- 
fields. 

• p(x) = YiiP( x i) : a (normalized) multidimensional 
probability density function, which is the product of the 
one-dimensional probability density functions p(x i ). 



• B v (x) : a product of the exponential one-body op- 
erators arising from the auxiliary-field transformation. 
From Eq. f), B v (x) = J], e x ^^ . 

• B(x) : the product of B v (x) with all other one-body 
exponential operators that do not depend on the aux- 
iliary fields x, and all the necessary scalar prefac- 
tors. For the projector in Eq. (|8}, B(x) = e ArET ■ 

With these notations, Vgs takes a generic form of a high- 
dimensional integral operator: 

V gs w / dxp{x)B(x) . (9) 
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D. Wave function representation 

We write our wave functions in terms of the basis functions 
\Xi)- A single -particle wave function is written as 

i^=E^i^>=E^ c iio)=^io). (10) 

i i 

A single-permanent, TV-Bosons wave function is given by 

\4>) =#$...^10). (ID 

In general, the exact ground state wave function is a superpo- 
sition of such permanents. Unlike the fermionic case, where 
the particles occupy mutually orthogonal orbitals, there is no 
such restriction on the orbitals here. We use this freedom in 
our method to have all the bosons occupy the same orbital in 
\4>), which greatly simplifies the computation Il4j| . We will re- 
fer to this as identical orbital representation (IOR). The most 



where 

V({x m , Vn}) = Ilm dXmWn d&n , 

P({x m ,y n }) = H m p(x™)I\ n p(yn) , 

and in the last line we have introduced the shorthand 

(v({xn})\ = (^ T \U n B(S n ); 
\my m }))^n m B(ym)\^ T )- 

The Metropolis simulation is carried out by sampling the 
probability density function defined by the integrand in the 
denominator. Given the choice of *f? T in the identical-orbital 
representation, this readily applies to bosons, which is how the 
model calculation by Sugiyama and Koonin Il4ll was done. 
The total length of the imaginary time is predetermined by 
At and the number of B operators in the product. 



III. NEW METHOD FOR BOSONS 

In this paper we formulate a new approach for ground-state 
calculations of bosons with branching random walks. There 



important virtue of this representation is that the exponen- 
tial of a one-body operator A transform a single-permanent 
wave function \<f>) into another single -permanent wave func- 
tion m 

e k \<l>) = W)> (12) 

In particular, B(x) in Eq. (I12t transforms a single permanent 
\(f>) into another single permanent \<fi'). (In Appendix lAl we 
include a brief summary of properties of wave functions in 
IOR.) 



E. Metropolis AFQMC 

Standard AFQMC calculations ffl employ Metropolis 
Monte Carlo to compute various ground-state observables, 



(13) 



are several advantages in implementing the Monte Carlo sam- 
pling as a random walk process. It is a true ground-state for- 
malisms with open-ended random walks which allow projec- 
tion to long enough imaginary-times. The sampling process 
can be made much more efficient than in standard AFQMC, 
by virtue of importance sampling with ^ T to guide the ran- 
dom walks. It also leads to a universal approach for bosons 
and fermions, where it is necessary to use the random walk 
formalism in order to implement a constraint to deal with the 
sign and complex-phase problems 

A key observation is that we can choose an IOR single- 
permanent wave function as the initial wave function |\I/ T ). 
At each imaginary timestep r = n At in the projection in 
Eq. Q, the wave function is stochastically sampled by a col- 
lection of single-permanent wave functions {|</>[ )}, where 
the index i (in Cursive letter) is different from the basis in- 
dex i. From Eqs. and t\2i . we see that, with each walker 
l^ "*) initialized to l^-p) in IOR, the resulting projection will 
lead to a superposition of single-permanent wave functions, 
all of which are in IOR. 

Each permanent evolves by the stochastic application of 
■Pgs, as follows: we randomly sample x from the probability 



■^gs|* T ) 



fD({x m ,y n }) P({S m ,Vn}) (grlTL B{x m ) A n» B(y n )\V T ) 
JV({x m , y n }) P({x m ,y n }) (* T |]l m B{x m ) Y[ n B{y n )\^ T ) 

(*/({*«})! iW{&})> 



$V{{x m ,y n })P({x m ,y n }) (r}({x m })\(j){{y n })) 



mSmWdyn})) 



JV({x m ,y n }) P({x m ,y n }) (v({xm})\<t>({yn})) 

I 
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density function p( x ), then apply B(x) on \ <p\ 



\4>\ 



(t+At) 



(14) 



We will call these permanents random walkers. The collection 
of these random walkers at each imaginary-time step is also 
referred to as population. 

The population must first be equilibrated so that the ground- 
state distribution is reached. After equilibrium the ground 
state is given stochastically by the collection of single- 
permanent wave functions 



l*o)=X>> 



(15) 



Measurement of ground-state observables can then be carried 
out. 

The random walk process naturally causes the walker's or- 
bitals to fluctuate. In order to increase sampling efficiency, 
we may associate a weight factor w\ to each walker \cf>i). For 
example, we can use the walker's amplitude as the weight fac- 
tor: 



A better definition of the weight will be introduced later when 
we discuss importance sampling. We duplicate a walker when 
its weight exceeds a preset threshold. Conversely, walkers 
with small weight (lower than a predetermined limit) should 
be removed with the corresponding probability. In this way, 
the walkers will have roughly the same weight. This results in 
a branching random walk. 



A. Measurement: "brute force" and mixed estimators 



(outside of the statistical errorbar) and converges to the exact 
value as iV^r is increased. In practice, however, the use- 
fulness of the 'brute force" estimator is limited to smaller 
systems. In general it will have large variances. Reducing 
the variance is expensive because (A) hf scales as C(iV^ lkr ), 
where N^k? is the size of the population used to represent 
l*o). 

The simplest approach to measuring the observables is the 
mixed estimator, i.e. 



(A 



(iM*o) 



(18) 



For example, to compute the ground-state energy, we can in- 
troduce the so-called local energy E h [ip T , <fi]: 



(19) 



The ground state energy is obtained from the weighted sum of 
the local energies associated with each walker: 



Ei( 7 M ( M £ L[y ; T>'^] 



(20) 



The local energy for each walker can be computed using the 
formula given in AppendixlAl 

The mixed estimator in Eq. i ll 8t is exact only if the operator 
A commutes with the Hamiltonian. Otherwise, a systematic 
error arises. Nonetheless the mixed estimator often gives an 
improvement over the purely variational estimator: 



(Ay 



(*tI4*t) 
(*tI*t) 



(21) 



Two formulas are often employed to correct for the systematic 



The ground-state value of an observable A is its expectation 
value with the ground-state wave function: 



(^)g,. 



($o|i|$o) 
($ol$o) 



(16) 



In principle, we can use the same Monte Carlo samples as both 
(<£>o| and |$o)- A "brute force" measurement on population 



{\4>\ )} at imaginary-time r is then given by 



\ A /b{ 



(17) 



and the estimator (A) hf is the average of such measurements. 
The "brute force" estimator is not useful in real-space based 
QMC methods such as diffusion Monte Carlo, because the 
overlaps between different walkers would lead to (5-functions. 
Here the walkers are non-orthogonal mean-field wave func- 
tions, and Eq. dl7> is well defined in principle. The estima- 
tor is exact for all observables in the limit of large -ZV w ik r . 
The ground-state energy estimated in this way is variational, 
namely, the computed energy lies higher than the exact value 



{A) atapl =2(A) mix -(A) T ; 



(A) extrap2 



(A) T 



(22) 
(23) 



The second formula is useful for quantities such as density 
profile, where it must be nonnegative everywhere. These cor- 
rections are good only if |^ T ) does not differ significantly 
from |$o)- In general, we need the back-propagation scheme 
to recover the correct ground-state properties. We will de- 
scribe this method after introducing importance sampling. 



B. Importance sampling 

In practice, the efficiency of the bare random walk de- 
scribed earlier is very low, because the random walks "ran- 
domly" sample the Hilbert space, and the weights of the walk- 
ers fluctuate greatly. This results in large statistical noise. We 
formulate an importance sampling procedure fl^.E2il — using 
the information provided by the trial wave function \^f T ) — 
to guide the random walk into the region where the expected 
contribution to the wave function is large. 
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I. Importance-sampled random walkers 



that the walkers propagate in the following manner: 



An importance-sampled walker also consists of a perma- 
nent and a weight, although the weight will be redefined ac- 
cording to the projected overlap of the permanent with the 
trial wave function. The purpose is to define a random walk 
process which will lead to a stochastic representation of the 
ground-state wave function in the form 



1*0) = E 



\4>i) 



(24) 



where Wi is the new weight of the walker. The overlap enters 
to redefine the weight factor such that walkers which have 
large overlap with |^ T ) will be considered "important" and 
will tend to be sampled more. Such walkers will also have 
greater contributions in the measured observables. Since the 
permanent now appears as a ratio \cf>i) / (\& T |(^j.}, its normal- 
ization is no longer relevant and can be discarded, unlike in 
the unguided random walk. The only meaningful information 
in \(f>i) is its position in the permanent space. 



2. Modified auxiliary-field transformation 

Now we describe the random walk process for the modified 
walkers. The goal is to modify "P gs in Eq. such that the 
random walk process leads to random walkers with the char- 
acteristics described above in Eq. J24> . The basic idea is the 
same as that in Ref. Il9| The main difference is that here we 
are dealing with bosons. In addition the HS fields in Ref. [3 
are discrete Ising-like, which allowed simplifications in the 
importance sampling, while here the auxiliary fields are con- 
tinuous and thus a more general formalism will be developed. 
Our mathematical derivation here follows that of Ref. |22l Up 
to now we have assumed that (^ T \(j)i) is real and positive. 
There is therefore no additional subtlety with the meaning of 
importance sampling and the correct form of the overlap to 
use, which Ref. |22| addressed in the context of fermionic cal- 
culations with general interactions. 

To derive the importance-sampled propagator, we plug 
Eq. (I24> into Eq. ( I3bt . We will focus on the two-body prop- 
agator, which is evaluated stochastically and is therefore af- 
fected by importance sampling in a non-trivial way. 

The modified propagator, V ga , consists of two parts. The 
first part is the transformation introduced in Eq. l|5}, which we 
now rewrite in the following form: 



1 



dx i 



/At {x—x)v 



(25) 



where we have added an arbitrary shift x to the auxiliary field 
x in the auxiliary-field operator. This is a change of variable 
in the integral on the right-hand side and does not alter the 
result of the integral. The new propagator 7^ must preserve 
the representation of | $o) m the form of Eq. (I24> : this dictates 



w 



(t+At 



m 



t+At) 



,(t) 



\<t>\ 



(t+At)> 



(26) 



From this requirement comes the second part of 
the modified propagator, which is the overlap ratio 



,(t+At), 



/<* t I# 



This factor is obtained by 



bringing the term (\l/ T |</>[ r+AT ') in Eq. i26\ to the right-hand 
side. It depends on an d the specific path in auxiliary- 
field space, and will "guide" the random-walk toward the 
region where (^f T \<pi) is large. 

Combining the two parts gives an importance-sampled 
propagator of the form 



where 



dxp(x)W(x, 4>)B{x — x) , 



(27) 



{^ T \B{x- x)\<t>) cS . s _i s . s (2g) 
(* T |0) 



is the aggregate of all the scalar prefactors in the modi- 
fied propagator. This propagator takes {w[ T \ {(f)^)} and 



advances the population to {w\ 



(t+At) I , (t+At) 



)}, both of 



which represent |$q) m the form of Eq. ( I24t . 

Monte Carlo sampling of the new propagator Vg S is similar 
to the one without importance samping. We sample x from a 
normal Gaussian distribution, and apply the operator B (x—x) 
to the current walker |</>[ r ). But now we accumulate an extra 

multiplicative weight factor W{£, 4>i) every time we apply 
Eq. CZ7): 



w[ t+At) ^W(x,^) w { T K 



(29a) 
(29b) 



Here we use the customary notation of vector dot product, e.g. 
x ■ x = Yl,ilk x i- Note that the weight factor W(x, (j)^) de- 

,(t+At)\ 



pends on both the current i 
positions. 



and future 



j ) walker 



3. The optimal choice for auxiliary -field shift x 

The optimal importance sampling is achieved when each 
random walker contributes equally to the estimator. We there- 
fore choose x to minimize the fluctuation in the weight factor 
Wi. The fluctuation in Wi will be minimized if we minimize 
the fluctuation in the prefactor Eq. d28l >. We do so by requiring 
the partial derivatives of this prefactor to vanish with respect 
to Xi at its average (Xi = 0): 



d 

dx. 



(tt T |S(£-s)|&) 
<* T |tH> 



0. 



x i= 
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It is sufficient to expand the exponentials in terms of At and 
require the term linear in Xi to vanish, since this is the leading 
term, containing \J At. The others contain higher-order terms 
and are vanishingly small as At — * 0. The best choice for x_ t 
that satisfies this requirement is 



(30) 



This choice depends on the current walker position as well 
as |^x)> which is t° be expected, since the objective for the 
shift is to guide the random walk toward the region where 




is large. With x determined, the algorithm for the 



random walk, as given in Eq. ( I29> . is now completely speci- 
fied. 



4. Local energy approximation 

We can furthermore approximate the prefactor W(x, 4>) in 
Eq. (I28> to obtain a more elegant and compact expression. 
After rewriting the prefactor in the form of an exponential, 
expanding B(x—x) in terms of At, and ignoring terms higher 
than O(At) in the exponent, we obtain 



n 



C -2 



iAT(l-I?)(5?-D?) iA™ 2 



where 



(* T ltflfr) 



(31) 



(32) 



The product is over the basis index i, which should be dis- 
tinguished from the walker index i. The latter is held fixed 
here. The first exponential in Eq. (13 1 1 can be ignored by not- 
ing that the average value of x\ with respect to the Gaussian 
probability density function is unity. Setting xf — > 1, i.e., 



evaluating the exponential at the mean value 



is justi- 



fied because vf and vf do not change drastically within one 

timestep. We also note that £\ v f = — (* T j V'|<0i)/(* T j<Ai)> 
which is the mixed-estimator of the potential energy with re- 
spect to the walker \4>i). Combining this term with the similar 
contribution from the kinetic propagator, we obtain a simple, 
approximate expression for Eq. J2 81 : 



W(x,^ ] ) 



Ar(B T -E L [vI/ 1 



(33) 



where E L [^ T ,^] is the local energy of tfii as defined in 
Eq. H9i . Note that, contrary to Eq. (I28> . this form depends 
only on the current walker position and not the future, al- 
though in practice a symmetrized version can be used which 
replaces the local energy by the average of the two. For a 
good trial wave function, the local energy fluctuates less in the 
random walk. If the trial wave function is the exact ground- 
state wave function, the local energy becomes a constant and 
the weight fluctuation is altogether eliminated. This bears a 
close formal resemblance to the importance-sampled difus- 
sion Monte Carlo method. 

The algorithm resulting from Eq. d33l is an alternative to 
Eq. |28i . The two are identical and exact in the limit At — > 0, 
but can have different Trotter errors. 



C. Measurement: back propagation 

With importance sampling, the mixed estimator in Eq. Jl 8I > 
is given by: 



4- 1 



(34) 



f/'i 



For example, the ground-state energy is 

As mentioned earlier, the normalization of cf>i is irrelevant be- 
cause <j>i only appears in ratios in any formula that defines the 
algorithm: Eqs. <ZEJ, OS, 03, and Eq. @U- We can 
(and should) normalize the permanent as needed, and discard 
the resulting normalization factor. 

The mixed estimator is often inadequte for computing ob- 
servables whose operators do not commute with the Hamilto- 
nian. In some cases the error due to this noncommutation is 
unacceptable. For example, the condensate fraction in the at- 
tractive trapped Bose-Hubbard model is greater than 100% if 
the Green's function (c\cj) is estimated using the mixed esti- 
mator. Therefore we have to propagate the wave functions on 
both the right- and the left-hand side of the operator: 



(tt T |e- T bP^i|$ ) 
(* T |e-^|<I> ) 



(35) 



This estimator approaches the exact expectation value in 
Eq. d 1 6I > as T bp is increased. Z hang and co-workers proposed 
a back-propagation technique 1 19] that reuses the auxiliary- 
field "paths" from different segments of the simulation to ob- 
tain ($q P | ee (<I> T |e~ Tb p ff , while avoiding the A^ lkr scaling 
of a brute-force evaluation with two separate populations for 
(<3>o| and |3>o)- Here we give a more formal derivation and 
description of the technique, and implement it to bosons. 

At imaginary-time t, the population is { 1 4>[ T ' ) }, which rep- 
resents |$o) in the form of Eq. ( I24t . The propagator in the 
denominator can be viewed equivalently as operating on the 
left or the right. The latter view is precisely the "normal" 
importance-sampled random walk from r to the future time 
t' = t + T bp , which consists of n bp ee T bp /Ar steps. 
We first assume that there is no branching (birth/death of 
walkers), i.e., the weights are fully multiplied according to 
Eq. (128b . The random walk of each walker will generate a 
path in auxiliary-field space. For convenience we will denote 



the path-dependent operator B[x\ T ' — x i^y')} by Bl T> , and 
weig ht factor W^, 4^ ] ) by w[ t) . Further we will denote 

~ (t) 

the time-ordered product of B^ from imaginary-time r to 

- ( T '- T ) ( t ) 

t by B\ , and correspondingly the product of by 

(t' -r) 

W> . Each path defines a product 



(t') 



(36) 
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Collectively these products give a stochastic representation of 



Replacing the operator e~ Tb p H in the numerator and de- 
nominator of Eq. J35i with Eq. 1361 . and using the expression 
for |$ ) given by Eq. ( 1241 . we obtain 



3. the population {|?y| Tbp ^)} is then generated by back- 
propagation using Eq. j39> : 

4. this population is matched in a one-to-one manner to 
{|<^ T ' ) )}, weighted by the weight at the later time, 



(A) 



£t<* 



bp 



■w, 



{t';t)£(t':t) £ „„(t) u (t) 



w 



, and the estimator is formed. 



Aw{ T >\cj>Y 



Ei(* T l 



<* T i^ T,) > 



(t' :t ) ^( r ' :T )„„( r )u( r ) 



or 



(37) 



In the back-propagation the propagators are, as shown in 
Eq. fl39l >. idential to those in the forward direction, but in re- 
verse order in imaginary-time. As in the normal walk, the nor- 

malization of \ r q i bp ) does not enter in the estimator. Similar 
to the mixed estimator, this procedure can be repeated period- 
ically to improve statistics. Evidently this estimator is exact 
in the limit of large r bp . 

We have assumed that there is no branching within the in- 
terval r bp . In practice, a population control scheme is often 
used which causes birth/death of walkers. This does not affect 
the derivation above or the basic algorithm. The effect on the 
implementation is that a list of ancestry links must be kept for 
the forward steps, which indicates the parent of each walker 
at each step in the imaginary-time duration r bp . As a result of 
branching, two or more (ry|'s may share the same segment of 

the paths in their "past" and the same parent |<?!>| T ' ) ). The esti- 
mator remains exact for large r bp . Branching or weight fluctu- 
ation does have a more serious practical implication, however. 
As r bp is increased, more and more (77 1 's will be traced back 

Note that each of these 77's originates from the trial wave func- to the same parent \<j>^). Or equivalently, fewer and fewer 



Using the propagation relation in Eq. d29l . we can show that 

i.e., the denominator in Eq. d37t reduces to J2i W [ T This 
result is to be expected, and can also be seen by completing the 
n bp steps of the "normal" random walk we discussed above. 
With importance sampling, the Monte Carlo estimate of the 
denominator is simply given by the weights at time r'. 

To simplify the numerator we associate a back-propagated 



wave function with each walker 



D 



( T + T b B : T ) 



(39) 



tion |^ T ), and is propagated by applying the B's in reverse 
order, as implied by the Hermitian conjugation. We may then 
write Eq. d37i in the following form: 



{Tbp, \m 



(r) 



(A) 



bp 



^ i 1 ( T bp)|j 1 W\ 



(t') 



(40) 



permanents in the set {|</>{ r ^)} will contribute to the estimator. 
This results in a loss of efficiency or an increase in variance. 
Better importance sampling will help improve the situation, 
often greatly, by reducing fluctuations in weights, although 
the problem will always occur at large enough r bp . In our ap- 
plications to date we have rarely encountered the problem and 
find that the computed observables converge quite rapidly (see 
section|V]fbr illustrative results). 



The estimators in Eqs. \35\ and d40l parallel that of the stan- 
dard AFQMC estimator in Eq. dl 3I >. The |</>)'s and (77 1 "s have 
similar meanings. The only difference lies in how the paths 
are generated. Here an open-ended random walk is used to ad- 
vance an ensemble of paths from r to t', which result in fluc- 
tuating weights that represent the path distribution. In stan- 
dard AFQMC a fixed length path (corresponding to r bp + r cq , 
with T cq being the minimum time for equilibriation or, failing 
that, the maximum time that can be managed by the calcula- 
tion) is moved about by the Metropolis algorithm, which elim- 
inates branching by the acceptance/rejection step. In other 
words, the estimators in Eq. dl 31 and Eq. d40l > are the same 
except for the weights. 

Eq. d40l > defines an algorithm for obtaining the estimate of 
(A) h via the following steps: 

1. A population is recorded as {|0| T ''}}; 

2. as the random walk continues, the path history is kept 



IV. TRAPPED BOSON GAS : MODEL AND 
IMPLEMENTATIONS OF QMC AND GP METHODS 

In this section we discuss the model we use to describe 
a single-species, Bose atomic gas with pair-wise contact 
interaction, confined in a harmonic trap in one- or three- 
dimensions. We then describe the implementations of both 
our QMC method and the standard mean-field GP approach 
to study this model. Numerical results will be presented in the 
following section, Sec.lVl 



A. Model 

We use an effective potential characterized by the low en- 
ergy atom-atom scattering length, a s . The two-body interac- 
tion takes a simple form 



for a time interval r, 



bp' 



Tit \ 47ra ^%/ 

U(ri - r 2 ) = — — — d(ri - r 2 J 



(41) 



m 
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For this effective potential to be valid, several assumptions are 
made; for example, the dominant effect is from s-wave scat- 
tering, and \a s \ is much smaller than the average inter-particle 
spacing. For more details we refer the reader to Ref.yl In the 
alkali gases these conditions are in general well met, and the 
model potential can be expected to give quatitative informa- 



tion, although care must be taken to validate the conditions. 

We now derive the Bose-Hubbard model from the standard 
many-body Hamiltonian of the trapped boson problem in d- 
dimension. In the continuous, real space, the Hamiltonian is 
given by: 



H = K + V = [ d 3 rSHr) (- J-V5! + hmuZr 2 ) $(r) 

J V 2m J (42) 



1 4na s h 2 

2 m~ 



d 3 r x / d 3 r 2 ft{r 1 )ft(r 2 )5{r 1 - r 2 )^(r 2 )^( ri ) 



The first term is the one-body Hamiltonian K, which consists 
of the kinetic energy and the (external) confinement poten- 
tial. V is the interaction Hamiltonian, which is the sum of all 
the two-body potentials. The characteristic trap frequency is 
u>o, which is related to the so-called oscillator length scale by 
a ho = y/h/mu; . 

We introduce a real-space lattice, with a linear dimension of 



L, in a simulation cell of volume (2ri,) d . The lattice spacing is 
therefore <; = 2r& / L. Further we will consider only a spheri- 
cally symmetric trap here for simplicity. We truncate the sim- 
ulation cell accordingly and assume that the wave function is 
negligible outside the maximum sphere enclosed by the cell. 
(Generalization to inhomogeneous traps is straightforward.) 
The discretized Hamiltonian corresponding to Eq. (1421 is 



l-'l T, c h 

^ jenn(i) 



- 2dc 



s |2„t. 



cjc-clc- 



(43) 



where cj and c, are the usual creation and annihilation opera- terparticle spacing, but larger than the scattering length: 
tors at site i. The Hubbard parameters t, U, and k are related 

to the real, physical parameters as follows: \a s \ < (< p~ x l d . (45) 

With negative a s , the particles tend to "lump" together due to 
the gain in the interaction energy. This is a situation where 
we especially have to be aware of the validity of the effective 
potential. As mentioned we will do a consistency check at 
the end of the calculation to ensure that the occupancy of the 
lattice points are less than unity. 

where for simplicity we have set h = m = 1. The lattice co- B. Implementation of QMC 
ordinate^ is related to the real coordinate by = (L/2ri y )r i , 

and f is the lattice coordinate of the trap's center. Note that Implementation of our QMC method for this model is 

a s is the true scattering length only in three-dimensional sys- straightforward. The number of basis M is equal to the num- 

tems. Nonetheless we will retain the symbol a s in Eq. d44bl ber of lattice sites inside the truncated sphere of radius r b . The 

as a convenient measure of the interaction strength in any di- two-body term in Eq. l|43) is in the desired form of Eq. 

mension. With a negative U, the HS transformation in Eq. @ leads to 

In the discretized model our resolution is limited by the lat- M aux iliary fi elds, with one-body propagators in the form of 

tice spacing. This is consistent with the conditions of validity exp(i/ Arjf/JxjnJ, where h i = is the density opera- 

of the model interaction in Eq. d41i . as it in a sense "inte- tor. Our trial wave function |^ T ) is the Gross-Pitaevskii (GP) 

grates out" the short-range dynamics. In this model our lattice wave function ^Qp, which we describe in the next subsection, 

constant ( must be much smaller compared to the average in- We mention here a technical point in the implementation. 



U 



1 

47ra s 



(44a) 
(44b) 
(44c) 
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The ground-state projection in our method involves the appli- 
cation of one-body propagator in the form of e A on a single- 
permanent wave function \<j>). This usually translates into a 
matrix- vector multiplication in the computer program, which 
generally costs 0(M 2 ). Often there are special properties of 
A that can be exploited to evaluate the one-body propagator 
more efficiently. In the Bose-Hubbard Hamiltonian, the only 
non-diagonal part of the Hamiltonian in real space is the ki- 
netic operator in K. We can separate it from the other one- 
body operators and apply the kinetic propagator in momentum 
space. Wave functions are quickly translated between these 
two representations using the Fast Fourier transform (FFT). 
In this way, the actual application of e~^ ArK involves only 
diagonal matrices; thus the overall cost for each e~^ ArK op- 
eration is reduced to 0(M log M). We observe in our calcu- 
lations that the additional Trotter error is much smaller than 
the error already introduced in the original breakup, Eq. (0J. 



the non-interacting Hamiltonian (with U = 0). The density is 
fed back to construct the initial Hamiltonian H G ° p in ( I48> . Di- 
rect diagonalization of this one-body Hamiltonian yields its 
ground state I'&Qp). We thus obtain an updated density n| 

and a better Hamiltonian -ff GP \ This procedure is iterated un- 
til the desired convergence criterion is satisfied. We choose 
our convergence condition to be: 



/rfr|y( t+1 >(r)-yW(r)| 
i/diV*+i)(r) + (^(*)(r)| 



(50) 



where e is a small number (usually on the order of 10~ 13 for 
double precision numbers). 

The second method we use to solve Eq. d48l > avoids the 
diagonalization procedure. It is closely related to the QMC 
method, both computationally and formally (see Sec. lVM . We 
use the ground-state projector e~ ATH <^ F : 



C. Implementation of Gross-Pitaevskii self-consistent equation 

The Gross-Pitaevskii (GP) wave function <j> GP is the single- 
permanent wave function 

$ GP (ri,r 2 , . . .rjv) = ¥>(ri)<p(r 2 ) ■ •■<p(r N ) , (46) 

which minimizes the expectation value of the ground-state en- 
ergy. Such a wave function satisfies the self-consistent Gross- 
Pitaevskii equation OS EH 



ft 2 

TV - 1 Aira.h 2 



VV(r) + ±mw 2 |r-r |V(r) 



(47) 



TV m 



\(p(r)\ 2 ip{r) = Hif(r). 



[We keep the prefactor (TV — 1)/TV, since we will study both 
large and small values of TV.] 

To compare our QMC results to those of mean-field, we 
carry out GP calculations on the same lattice systems. The 
discretized GP Hamiltonian in the second-quantized form is: 



E ( E 4c, - 2d4c) 



if. -f n i 2 4c. 



l K E ,A 

i 

TV - 1 



(48) 



TV 



^E 



- t 1-2 N 



Here n, is the expectation value of the density operator: 

J, 



($GPK C »I $ GP) 
(*GPI*GP) 



(49) 



We have implemented two methods for solving the GP 
equation. The first is the usual self-consistent iterative ap- 
proach. We generate an initial density profile, n\°\ by solving 



(e 



-ArH r _ 



(0)\ 



1$ 



GP/ 



(51) 



The initial wave function is arbitrary and can be, for example, 
chosen again as the solution with [7 = 0. The feedback mech- 
anism through the density profile h i remains the same. By 
using the same Fast Fourier transform for the kinetic propaga- 
tor as described in subsection llVBI a speed gain is obtained, 
especially for large systems. In practice we have often found 
this method to be a simpler and faster alternative to the first 
method of diagonalization and iteration. Note that the scalar 
term — | U , does not affect the projection pro- 
cess, but with it H GP corresponds to the original many-body 
Hamiltonian in that ($ GP |7J GP |<I> G p) = ($ GP |H|$ GP ). 



V. RESULTS 

In this section we present results from our QMC and GP 
calculations in one-, two-, and three-dimensions. To validate 
our new QMC method and illustrate its behavior, the major- 
ity of the calculations will be on systems where exact results 
are available for benchmark. These include small lattices, 
which can be diagonalized exactly, and the case of attractive 
(5-function interactions in one dimension, where analytic solu- 
tions exist. For the purpose of presenting the method to facili- 
tate implementation, some numerical results and comparisons 
are shown in detail to illustrate the behavior and characteris- 
tics of the method. 

Most of the results we present here will be for attractive 
interactions, where the method is exact and is free of any 
phase problem ll22ll from complex propagators (see subsec- 
tion lV C> . Such systems therefore provide a clean testground 
for our new method. In addition, with attractive interactions 
the condensate in 3-D is believed to collapse beyond a critical 
interaction strength or number of particles. Mean-field cal- 
culations l30ll estimate the collapse critical point to be about 
Na s /a ho = —0.575. The exact behavior of the condensate 
near the critical point is, however, not completely clear, as 
many-body effects are expected to have an impact. At the end 
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of this section we will also show some preliminary results for 
larger systems with both attractive and repulsive interactions 
in 3-D. 

We measure the ground-state expectation values of the fol- 
lowing quantities: the ground-state energy, kinetic energy 
(T), external confining potential (Vtrap), interaction energy 
(^2b)> density profile (n i ), and the condensate fraction (of- 
ten abbreviated "cond.frac." in the tables and figures). The 
condensate fraction is defined as the largest eigenvalue of the 
diagonalized density matrix J3|. If we write the one-body 
Green's function matrix (c\cj) in terms of its eigenvalues 
{n a } and eigenvectors {x Q («)}: 



then the largest eigenvalue divided by the total number of par- 
ticles gives the condensate fraction. 



to rather large At values. We see that both quantities converge 
to the exact results as At — > 0. 

Trap Potential Energy 



0.850 



0.840 
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0.820 



0.810 



0.800 



QMC (T BP = 4.0) 
Exact diag. 




A. Comparison with exact diagonalization: a s < 

The many-body Hamiltonian (I43> can be diagonalized ex- 
actly for small systems to benchmark our QMC calculation. 
We compare our QMC results with exact diagonalization for 
a one-dimensional lattice of 13 sites, and study its behavior 
for different values of the interaction strength a s and number 
of particles N. 

The first system we study has 5 bosons, with t = 2.676, 
U = —1.538, K = 0.3503. These values correspond to the 
physical parameters a ho = 8546 A and a s = —5.292 x 
10~ 6 A~ x . (Recall that, by our definition, a s in 1-D does 
not have the dimension of length, and is not the scattering 
length itself.) Table H] shows the comparison of the quantities 
computed using three methods: QMC, GP, and exact diag- 
onalization (ED). The statistical uncertainty of QMC results 
are presented in parantheses. We see that the agreement be- 
tween QMC and ED is excellent. GP makes significant errors 
here because of the sizable interaction strength as well as the 
small number of particles. 

TABLE I: Comparison of QMC calculation against exact diagonal- 
ization (ED) and Gross-Pitaveskii (GP). The system has 13 sites, 5 
particles, t = 2.676, U = -1.538, K = 0.3503. In the QMC cal- 
culation we use Ar = 0.01, r bp = 4.0, and the GP solution as the 
trial wave function. 

Type g.s. energy (T) (Vtrap) (Vzb) cond.frac. 

ED -1.009 4.278 0.8427 -6.129 95.59% 

QMC -1.008(2) 4.279(3) 0.8423(5) -6.129(2) 95.59% 

GP -0.493 3.919 0.7504 -5.162 100% 



To illustrate the convergence in imaginary-timestep At, we 
show in Fig. ^ the total energy and the average trap energy 
(Vtrap)- The former can be obtained exactly from the mixed 
estimator while the latter requires back propagation. To show 
the Trotter error, we have deliberately done the calculations up 



FIG. 1: Convergence of QMC observables with At. The system has 
the same parameters as in TableQ Exact results are shown as dotted 
lines. Lines connecting QMC data are to aid the eye. 

To illustrate the convergence of observables in back- 
propagation length, we show in Fig.|2]the various observables 
computed by QMC as a function of T bp . Separate calculations 
were done for different values of T bp . For all calculations, 
a small At value of 0.01 was used. We see that all quanti- 
ties converge to the exact results rather quickly, by T bp ~ 2. 
(The total energy (H) is of course exact for any T bp , including 
T bp = 0.) As we see from the energy expectations, this is in 
fact a system with significant interaction effects. Alkali sys- 
tems at the experimental parameters often have significantly 
weaker interaction strengths and the convergence rate is ex- 
pected to be even faster. 

Our QMC method is exact and therefore independent of the 
trial wave function ^> T , except for convergence rate and sta- 
tistical errors. In Fig. [3] we show QMC results obtained us- 
ing two different ^ T 's, the noninteracting solution and the 
GP wave function. The convergence of condensate fraction 
and trap energy are shown versus back-propagation time T bp 
for a system of 6 particles on 13 sites. The calculations lead 
to the same results. The quality of ^> T , however, does af- 
fect the variances of the observables and their convergence 
rates with T bp . For example, the noninteracting wave func- 
tion, which disregards the two-body interaction, is more ex- 
tended (in its density profile) than GP. Its mixed estimator is 
therefore worse than that with the GP trial wave function. The 
mixed-estimator for the ground-state energy is exact in both, 
but the variance is slightly larger with the former. 

We now show results for different systems with N from 2 to 
9 bosons, and varying interaction strengths. We note that if we 
keep the product U x (N — 1) constant, the Gross-Pitaevskii 
equation predicts the same per-particle energies and densities. 
For brevity, we shall refer to the curve in which U x {N — 1) 
is constant as the GP isoline. Deviation from the GP isoline 
is therefore an indication of the effect of many-body corre- 
lations. In order to show results on multiple systems at the 
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FIG. 2: Convergence of the computed observables versus r bp . The 
system is the same as in Table The different panels show five 
different observables. The horizontal axes are the back-propagation 
length. Exact results are shown as dotted lines, while GP results as 
dash-dotted lines. Solid lines are present only to aid the eye. 
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FIG. 3: Independence of QMC results on trial wave functions ("GP" 
for Gross-Pitaevskii, "nonint" for noninteracting solution). The sys- 
tem is the same as in Tabled except that here we use 6 particles. The 
horizontal axes are the back-propagation length. Lines connecting 
QMC data points are present only to aid the eye. 



same time we will scan GP isolines. Figure[4]shows the QMC 
and GP results as a function of the number of particles. In the 
GP calculations the per-particle quantities are constants. The 
QMC results, on the other hand, capture the effect of corre- 
lation. Both the total energy and the interaction energy are 
lowered from the GP results. The exact results deviate from 
GP more as the system becomes more correlated along the GP 
isoline, i.e. when U is increased or when N is decreased. Al- 
though N is too small here because of the limitation of ED, 
the results are representative of the general trend in larger sys- 
tems (see below). 

Figure [5] further illustrates the effect of particle correla- 
tion in this system. Although the exact interaction energy is 
lower than that of GP, the exact density profile is more ex- 
tended. This is also manifested in the average trap potential 
energy (Vtrap) /N, where the QMC results are 0.1981(8) and 
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FIG. 4: Comparison of QMC, GP, and ED results for different sys- 
tems. Calculations were done along a GP isoline U x (N — 1) = 
— 2.30t for up to nine particles in 13 sites. The graphs show the to- 
tal and interaction energies per particle. QMC and exact results are 
indistinguishable. GP is accurate in the limit of weak correlation but 
deviates more from the exact results as the system becomes more 
correlated. The solid lines are to aid the eye. 



0.1605(2) for N = 2 and 9 particles, respectively, while the 
GP value is 0.1501. In GP, interaction energy is lowered by 
increasing particle overlap, namely by shrinking the profile. 
In reality, the particles find a way to reduce interaction with- 
out statically confining to the central sites, resulting in a more 
extended one-body profile. 




5 6 7 

Lattice site 

FIG. 5: The normalized density profiles as an illustration of parti- 
cle correlation effects. Results are for 13-site systems along the GP 
isoline U x (N — 1) = — 2.30t. The normalized GP curve is iden- 
tical for any number of particles along this line. QMC results are 
shown for N = 2 and N = 9. The QMC results have very small 
errorbars and are indistinguishable from ED (not shown). The QMC 
density profiles are more extended, although the interaction energies 
are lower than GP, as shown in Fig. [4] 



B. Comparison with analytic results in 1-D: a s < 

The problem of an arbitrary number of untrapped bosons 
interacting with an attractive (^-potential in one dimension can 
be solved analytically l3lll . yielding analytic expressions for 
the total energy and density profile. In this section we carry 
out QMC and GP calculations and compare our results against 
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these analytic results, on systems of up to 400 bosons. The 
Hamiltonian in the continuous real space is 



N 



N 



H 



i ^ d 2 l ^ 



Xj ). 



(52) 



i>j=i 



The interaction constant (g > 0) is related to our Hubbard pa- 
rameters by g = \U / \ft |. The ground state of this Hamilto- 
nian is an iV-boson bound state. By fixing the center of mass 
at x — 0, we can eliminate the contribution from its overall 
motion, which leads to the following analytic expressions for 
the density profile ll32ll . 



N-l 
n=l 

and the total energy, 



n(N\) 2 e~ gnN ^ /2 
(N + n- l)\(N-n- 



E = -±g 2 N(N 2 - 1) 



(53) 



(54) 



In our QMC calculations, we again put the system on a real- 
space lattice. The lattice size is chosen to be large enough so 
that discretization errors are comparable to or smaller than sta- 
tistical errors. As the ground state of the system is a droplet 
in the absence of the external confining potential, the center 
of mass can slide in the calculation due to random noise. We 
therefore need to subtract the center-of-mass motion. Tech- 
nically, this can be accomplished conveniently in the random 
walk by treating the system with respect to its center of mass. 
In Appendix lEl we describe our method for this correction, 
which is applicable in any situation where the center of mass 
and relative motions need to be separated. In our calculations, 
the correction affects the kinetic and total energies as well as 
the density profiles. The results shown below were all ob- 
tained with such a correction applied. 

We first study a system of 20 particles with g = 0.154. 
Table HT1 shows the energies, and Fig. [6] the density profiles. 
This is a system where mean-field makes significant errors. 
Our QMC results are in excellent agreement with the exact 
results. 



TABLE II: Comparison of QMC and GP results to available exact 
results. The system has 20 particles and g = 0.154. A lattice of 1024 
sites was used, with At = 0.01 and r bp = 2.5. 



Type 



g.s.energy (T) 



(Van) 



cond.frac. 



Analytic result 

QMC 

GP 



-1.971 

-1.964(8) 2.044(8) 
-1.784 1.776 



-4.007(4) 99.76% 
-3.561 100% 
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FIG. 6: Comparison of calculated density profiles from QMC and 
GP with analytical results. The densities are normalized. The QMC 
errorbars are displayed every five data points to avoid cluttering the 
plot. The QMC profile is given by the dotted curve. The inset shows 
the same curves with logarithmic vertical scale, indicating that at 
large distances the density is exponential. 



as N is decreased, mean-field results deviate more and more 
from the exact results. For example, as we go from g = 0.01 
(N = 400) to 10 times the strength along the isoline, the sys- 
tematic error in the GP total energy increases roughly from 
0.5% to 5%. 
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FIG. 7: Comparison of the energy from QMC (crosses) with the 
exact answer (dotted curve) for different number of particles. Energy 
per particle is shown along the GP isoline g x (N — 1) = 4.0. The 
GP result is the flat, dash-dotted line. We use a lattice of 1024 sites, 
Ar = 0.01 and r bp = 4.0. 



We next scan systems with various numbers of particles by 
following the GP isoline g x (N — 1) = 4.0. The energy 
per particle is shown as a function of N in Fig. [7] for up to 
400 particles. Fig.|8]shows the density profiles for up to 100 
particles. Again, the agreement between QMC and exact re- 
sults is excellent. As the interaction strength g is increased or 



We now study the system along a different line, holding the 
interaction strength g fixed while scanning the number of par- 
ticles, again up to N = 400 particles. Figure [9] shows the 
behavior of (H)/N 3 for up to 400 particles, with g = 0.0403. 
At large N, the total energy is roughly proportional to A^ 3 . 
Compared to Figs. and [8] the interaction strength here is 
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FIG. 8: Comparison of the density profiles from QMC and GP with 
analytic results. The normalized densities are shown along the GP 
isoline g x (TV — 1) = 4.0 for several TV values. The system is the 
same as that in Fig.0 The GP density is the same for any TV on the 
isoline, and is given by the dash-dotted line. 



stronger at larger TV and weaker at lower TV, with the crossover 
at TV ~ 100. Most of the calculations are therefore more 
challenging numerically. Again QMC was able to completely 
recover the correlation energy missed by GP. At large TV, 
smaller timesteps were used and more computing was nec- 
essary to reduce the statistical errors. (Note that the errorbars 
appear larger at smaller TV in the plot because of the division 
by TV 3 .) 
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FIG. 9: Comparison of computed ground-state energy for different 
numbers of particles TV. The interaction strength is held constant at 
g — —0.0403. The total energy divided by TV' ! is shown as a function 
of TV for QMC, GP and exact calculations. Conservative parameters 
were used, with r bp = 4.0 in all case, and At = 0.01 for TV < 200 
and At = 0.005 otherwise. 



with 13 sites and 4 particles. The agreement between QMC 
and exact result is excellent. Results from GP are also shown. 
The GP and QMC density profiles have roughly the same size, 
as evident from the values of (Vtrap)- However, GP overesti- 
mates the interaction energy because it does not take into ac- 
count the particle-particle correlation. In the mean field pic- 
ture, expanding the density profile is the only way to lower 
the interaction energy, so that the particles overlap less with 
each other. (Note that (Vtrap) is indeed slightly larger for 
GP.) In reality, particles can avoid each other more effectively 
by means of many-body correlation. The QMC correctly re- 
covers this correlation, which lowers the total energy without 
spreading the density as much as GP does. 



C. Comparison with exact diagonalization: a s > 

We have shown that our new QMC algorithm is exact and 
works well for a wide range of systems with attractive inter- 
actions. If the interaction is repulsive (a s > 0, or equivalently 
U > 0) the one-body propagators resulting from t he HS trans- 
formation become complex, in the form of exp(iV AtU Xihj). 
The same algorithm applies in this case as well. In principle 
the complex one-body operator only requires a change to the 
corresponding complex operations. But in practice a serious 
phase problem occurs, which causes the calculation to lose ef- 
ficiency rapidly at larger interaction strengths. We discuss this 
problem and how to control it below. Our initial studies indi- 
cate that, for moderate interaction strengths, the algorithm as 
is remains very efficient and gives accurate results, allowing 
reliable calculations for parameters corresponding to experi- 
mental situations in 3-D. 

We benchmark our algorithm in one- and two-dimensional 
systems with repulsive interactions against exact diagonaliza- 
tion. Table [H]] shows results for a one-dimensional system, 



TABLE III: Comparison of QMC results against exact diagonaliza- 
tion (ED) and Gross-Pitaveskii (GP) in 1-D. Here we use 13 sites and 
4 particles; t = 2.676, U = +1.538, k = 0.3503; At = 0.01 and 
r hp = 2.5. 

Type g.s.energy (f) (Vtrap) (Vjffl) cond.frac. 

ED 4.24 1.18 1.793 1.269 98.5% 

QMC 4.24(2) 1.18(2) 1.790(8) 1.273(8) 98.6% 

GP 4.43 1.03 1.800 1.599 100% 



Table IIVI shows results for bosons in a two-dimensional 
trap, using a 4 x 4 lattice. The GP solution also exhibits the 
same behavior as in the 1-D calculation, in that the density 
profile is slightly more extended, and the interaction energy is 
overestimated. As in other cases, the QMC statistical errorbar 
on the condensate fraction was not computed directly, but we 
estimate it to be on the last digit. 

As mentioned earlier, the only modification necessary to 
the algorithm in order to treat repulsive interactions (a s > 0) 
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TABLE IV: Comparison of QMC calculations against exact diag- 
onalization (ED) and Gross-Pitaveskii (GP) projection in a 4 x 4 
lattice, with 4 bosons, t = 0.2534, U = +0.3184, k = 3.700; 
At = 0.01 and r bp = 2.5. 

Type g.s.energy (f) (Vtrap) (V2b) cond.frac. 

ED 6.000 1.818 3.8326 0.350 97.8% 

QMC 6.005(6) 1.817(2) 3.8325(2) 0.355(5) 97.8% 

GP 6.067 1.763 3.8359 0.469 100% 



exact ground-state wave function. We see that this is indeed 
the case in Table IV1 The interaction energy is lowered in the 
many-body calculation as expected. Interestingly, the external 
potential energy is lower than in GP. Consistent with this, the 
exact density profile is tighter than in GP, as shown in Fig. 1101 
The trend here appears different from what we observed in 
small 1-D trapped systems in Fig.|5J but consistent with the 
large untrapped systems in Fig. [8] We are presently carrying 
out more calculations to cover a wider range of parameters 
and study the role of dimensionality. 



is to allow complex arithmetic. A more serious problem can 
occur, however. The orbitals and the walker weights become 
complex numbers. Asymptotically the phase of these weights 
will be uniformly distributed in the complex plane. The de- 
nomitors in Eqs. (I34> and J40i will be dominated by noise, 
causing the Monte Carlo sampling efficiency to decay and ul- 
timately destroying the algebraic sca ling of QMC. This is the 
so-called sign or phase problem Jl9ll22ll . In real-space meth- 
ods this problem is connected to fermions, but here we have 
a situation where a phase problem appears in the ground state 
of a bosonic system. Physically, it is easy to see why a phase 
problem must occur. Our many-body wave function is being 
represented in IOR, with only one orbital in each walker. With 
a repulsive interaction, the only way to reflect correlation ef- 
fects, i.e., particles avoiding each other, is to make the orbitals 
complex. 

As we see below, our algorithm remains efficient and gives 
accurate results for large systems with scattering lengths cor- 
responding to experimental situations in 3-D. As the interac- 
tion strengths become much stronger, the phase problem will 
ultimately make the approach ineffective. We have done pre- 
liminary calculations in which we control the phase problem 
by applying a phaseless formalism described in Ref. |22| Our 
results indicate that the systematic errors introduced by the 
phaseless approximation are small for moderate interaction 
strengths. We expect to therefore be able to obtain accurate 
and reliable results for scattering lengths well into the exper- 
imental 'strong-interaction' regime achievable by Feshbach 
resonnance. 



D. Realistic calculations in three-dimensions 

In this section we present some test results on realistic 
systems of trapped particles in three-dimensions. QMC re- 
sults were obtained with back-propagation and conservative 
choices of At and convergence parameters. We expect the 
QMC results to be exact. We also carry out the corresponding 
Gross-Pitaevskii calculations, and make comparisons against 
our exact QMC results. 

TablelVlshows the result of a QMC calculation for 175 par- 
ticles in a three-dimensional trap. We choose a trap with a 
characteristic length a ho = 8546 A. The trap was discretized 
into a 15 x 15 x 15 lattice, in a range that corresponds to 
about 5 x a ho . The scattering length is a s = —22.4 A. In 
this regime the GP solution is a good approximation to the 



TABLE V: Comparisons of QMC and GP calculations for 175 parti- 
cles in a 3-D spherical trap, with a s — —22.4 A and a ho = 8546 A. 
The energies are displayed as per-particle quantities. Both the QMC 
and GP results are extrapolated to At — > 0. 

Type g.s.energy (T) (Vtrap) (Vm) cond.frac. 
QMC 16.979(6) 16.47(5) 6.54(1) -6.03(4) 99.73% 
GP 17.115 15.60 6.77 -5.25 100% 
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FIG. 10: Comparison of density profiles from the QMC and GP for 
175 particles. The system is the same as described in Table Ivl The 
QMC profile is more peaked and tighter than GP. 

We now turn to bosons with repulsive interactions in three- 
dimensional trap. We again use al5xl5xl5 lattice, and sim- 
ulate 100 bosons. We choose a scattering length a, of 80 A. 
This value is close to the experimental 39 K singlet 133] or 87 Rb 
triplet 1 34] scattering lengths. In Table IVT1 we show the cal- 
culated energies and condensate fraction. For this interaction 
strength, the impact of the phase problem on the statistical er- 
ror is small, and the QMC calculation is very efficient. The 
true condensate is, like in the 1-D repulsive case, tighter than 
that predicted by GP, with lower interaction energy. 

VI. DISCUSSIONS 

A. Connection between QMC and Gross-Pitaevskii projections 

The QMC method we have presented allows us to go be- 
yond mean-field and treat many-body effects. On the other 
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TABLE VI: QMC calculation of 100 particles in a three-dimensional 
trap. A lattice of 15 x 15 x 15 was used. The parameters correspond 
to a ho = 8546 A and a s = 80 A. The quantities displayed are for 
per particle. 



Type g.s.energy (f) 



(Vtrap) (V2b) cond.frac. 



QMC 24.687(9) 9.573(9) 11.933(5) 3.181(3) 99.80% 
GP 24.922 9.281 12.028 3.612 100% 



hand, it has a deep connection with the GP mean-field ap- 
proach. Our approach uses an HS transformation which leads 
to integrals of single-particle operators over auxiliary-fields. 
The mean-field solution can be regarded as the leading term 
in the stationary-phase asymptotic expansion of the exact so- 
lution l35ll . Our method evaluates this exact solution, which 
is in the form of many-dimensional integrals, by Monte Carlo. 
In this section we further comment on the formal connection 
between our importance-sampled QMC and the GP as done 
by projection (the second of the two GP methods discussed in 
subsection llV C> . 

Let us reconsider the two-body propagator in the modified 
AF transformation Eq. d25i . Let us suppose that we are now 
taking our first Monte Carlo step, where our walker is \<fi), and 
we will also use the same wave function as |\& T ). Following 
the discussion of the optimal choice of x in the same section, 
IIII Bl we know that x = is a stationary point with the choice 



(55) 



We can approximate the integral in Eq. d25i by the value of 
the integrand at $ = 0, which can be justified in the limit of 
small At. More explicitly, as At — ► 0, the Gaussian function 
becomes the most rapidly varying term in the integrand. To 
exhibit the asymptotic behavior of this integral, we change the 
integration variable to y = \J At x, so that the large parameter 
1/At appears in the Gaussian's exponent: 



e- 



Atv 1 



-At(|s 2 



dy 



r J/ 2 /2Ar 

V2ttAt 



,y(v-v) 



The dominant contribution to the integral comes from the 
maximum of the Gaussian function at y = 0. The asymptotic 
leading term of the importance-sampled many-body propaga- 
tor is therefore: 



(56) 



where K is the one-body term in the original Hamiltonian. 
Under this approximation, our random walk becomes deter- 
ministic, needing only one walker. If for the next step we use 
the updated wave function \<f)') to evaluate the new {-D.;} in 
Eq. d55> . we obtain a self-consistent projection with one-body 
propagators. In fact, the one-body Hamiltonian in the expo- 
nent of Eq. d56l is precisely the mean-field Hamiltonian. For 



example, for Bose-Hubbard model the last two terms in the 
exponent lead to the GP mean-field potential 



(57) 



Apart from the factor (N — 1)/N which approaches unity in 
the limit of large N, we have recovered the GP propagator. 
The projection with Eq. d56t lowers the variational energy for 
any initial \<fi) and is stationary when \<p) is the GP solution. 
This is why GP is the best variational wave function that has 
the form of a single permanent, and hence a reasonable trial 
wave function to use for most of our QMC calculations. 

It is also clear from the discussion above that the impor- 
tance sampling formalism allows us to have an optimal form 
of HS transformation, in that the HS propagator e v ^ v ~ vS> in- 
volves only the difference v — v. In other words, although 
in Eq. |7} we write the decomposition for the bare inter- 
action term, the importance sampling transformation effec- 
tively introduces a mean-field background based on the trial 
wave function and allows the HS to deal with only a residual 
quadratic interaction term, (v — v) 2 . 

To summarize, our QMC method reduces to GP if we eval- 
uate the many-body propagator by the stationary-point ap- 
proximation, using only the centroid of the Gaussian. The 
full method evaluates the many-dimensional integral over 
auxiliary-fields exactly by Monte Carlo. It captures the in- 
teraction and correlation effects with a stochastic, coherent 
ensemble of mean-field solutions. The structure of the cal- 
culation can be viewed as a superposition of the GP projec- 
tions that we have described. Our method therefore provides a 
way to systematically improve upon GP while using the same 
framework. 



B. Computing 

Because of the structure of QMC as a superposition of GP 
projections, our method scales gracefully with system size. 
As discussed in Sec. II V Bl the bulk of our method scales as 
0(M log M), with the significant speedup from using Fast 
Fourier transform. For example, the QMC calculation shown 
in Table lVII reauired less than 8 hours on a single Alpha EV67 
processor. The 1024-sites QMC calculation shown in Table ITU 
took about four hours to get good statistics, with very con- 
servative choices of At and other convergence parameters. 
It required about 1.3 gigabytes of memory, largely because 
of back-propagation path recording. In contrast, treated fully, 
the latter problem would mean the diagonalization of a sparse, 
Hermitian matrix containing (8 x 10 41 ) 2 elements. Although 
this can be reduced by exploiting symmetries, exact diagonal- 
ization of this problem is clearly not within reach with com- 
puting capabilities in the foreseeable future. 

We typically use hundreds of walkers in our calculation. 
The stochastic nature of QMC means the number of walkers 
fluctuates due to branching and killing of walkers with very 
large and very small weights (see subsection lllli. The popu- 
lation therefore must be controlled to ensure that it does not 
grow or decay too much, and that the walker weights have a 
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reasonable distribution. Our method to control the population 
is similar to that discussed in Ref. |3 

We comment on the effect of the number of particles, N, 
on scaling. Because of the use of IOR, the number of parti- 
cles does not enter in the propagation. It would then seem as 
though the algorithm might have a super-scaling in N. This 
is not true, of course, since the projector e~ ArH depends on 
N. For example, the shift vi has a factor of N in front (see 
Appendix[XJ, and the local energy scales with N. As a result, 
a smaller time-step must be used for larger N . The above ar- 
guement suggests a linear reduction in At as N is increased, 
which we have used as a rough guideline in our calculations 
to select the range of At to use. Extrapolations with separate 
calculations using different At values are then carried out. 



C. Conclusion and Outlook 

In conclusion, we have presented a new auxiliary-field 
QMC algorithm for obtaining the many-body ground state 
of bosonic systems. The method, which is based upon the 
field-theoretical framework and is essentially exact, provides 
a means to treat interactions more accurately in many-body 
systems. Our method shares the same framework with the GP 
approach, but captures interaction and correlation effects with 
a stochastic ensemble of mean-field solutions. We have il- 
lustrated our method in trapped and untrapped boson atomic 
gases in 1-, 2-, and 3-dimensions, using a real-space grid as 
single-particle basis which leads to a Bose-Hubbard model 
for these systems. We have demonstrated its ability to obtain 
exact ground-state properties. We have also carried out the 
GP mean-field calculations and compared the predictions with 
our exact QMC results. Our method is capable of handling 
large systems, thus providing the possibility to simulate sys- 
tem sizes relevant to experimental situations. We expect the 
method to complement GP and other approaches, and become 
a useful numerical and theoretical tool for studying trapped 
atomic bosons, especially with the growing ability to tune the 
interaction strengths experimentally and reach more strongly 
interacting regimes. 

From the methodological point of view, more work remains 
to be done with the repulsive case to deal with the phase prob- 
lem. We have shown that our method as it stands can be 
very useful for moderate interaction strengths. For stronger 
interactions, our preliminary study indicates that the phase- 
less approximation 12211 . which eliminates the phase problem 
but introduces a systematic error, is very accurate for scatter- 
ing lengths well into the Feshbach resonnance regime. We are 
currently examining this more systematically to quantify the 
extent of the bias. Because of the simplicity of these bosonic 
systems compared to electronic systems, they provide an ideal 
testbed, where for small sizes the problem is readily solved by 
exact diagonalization. 

A variety of applications are possible. The ground state 
of the Bose-Einstein condensates with both attractive and re- 
pulsive interatomic interactions can be studied for various in- 
teraction strengths, including the strongly interacting regime 
reached by Fesbach resonance. They can also be studied in 



different dimensions and under different conditions. In partic- 
ular, it would seem straightforward to generalize our present 
framework to study rotations and vortices, since we are al- 
ready dealing with complex propagators and wave functions 
in the repulsive case. In addition, it will be interesting to 
treat boson-fermion mixtures with our approach. As men- 
tioned, the auxiliary-field method is already widely used to 
treat strongly interacting fermion systems. 
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APPENDIX A: IDENTICAL-ORBITAL REPRESENTATION 

In this appendix we show that the matrix representation of 
an -boson wave function in AFQMC can be made particu- 
larly simple. In fermion calculations, we must use anMxJV 
matrix to represent a determinant, because the orbitals must be 
mutually orthogonal. In the boson case, however, this restric- 
tion is absent. The most general form of a many boson perma- 
nent is expensive to compute, having complexity of O(JVMI). 
But we can choose to make all the orbitals identical. In ma- 
trix language, we will have only an Af-row column vector. 
We will term this representation identical-orbital representa- 
tion — IOR. Each many-boson wave function in IOR has the 
form of a GP mean-field solution. Two conditions are neces- 
sary for this choice to be viable in the QMC: that an initial 
trial wave function of this form is allowed and that successive 
projections preserve the form. The only requirement for the 
former to hold is that the wave function in IOR not be orthog- 
onal to the true many-body ground state, and it is straightfor- 
ward to show that Eq. (1121 holds for a \<f>) in this form. More 
complex wave functions can always be generated by a linear 
combination of such wave functions. In fact, this is what we 
accomplish through our Monte Carlo simulation. 

In operator language, a single A -boson wave function \<j>) 
is given by 

|^)=^t..^t| 0) = ^t)^| 0)) 

N 

where <ft = ^Zq, c^q,. In matrix form, \<f>) would be M x N 
matrix cf> whose columns are identical. The overlap of two 
such wave functions is given by 

(V#) = Per (tp T ■ 4>) 
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where the bold-phased symbols if) and <j> represent the single- 
column vectors for ip and <f>, respectively. Similarly, for any 
one-body operator A, 

= N\ • A • 0)(*/> f • 0) Ar " 1 , (Al) 

where A is the matrix for A. The matrix element of a quartic 
(two-body) operator is given by: 

(Mbibib^M = n\ n(n - i)r a r^M^ ■ <t>) N - 2 . 

(A2) 



APPENDIX B: DROPLET CENTER-OF-MASS 
CORRECTION 

1. Correcting the density broadening 

To handle the droplet system given by the translationally 
invariant Hamiltonian in Eq. ( I52i . an extra ingredient is nec- 
essary in addition to the "basic" QMC algorithm that we have 
described. In a deterministic calculation, for example in GP, 
the motion of the center-of-mass (CM) can be simply elimi- 
nated by fixing it at the origin, as in Eq. J53t . In the QMC 
calculation, however, the orbitals fluctuate as they are propa- 
gated by B(x — x), where the random fields x are drawn from 
a Gaussian probability density. Random noise will inevitably 
cause the CM of the system to slide, undergoing a free diffu- 
sion whose average position is the origin. 

Left unchecked, this spurious CM motion will lead to an 
artificial broadening of the density profile. To correct for it 
in the density profile, we could simply shift the CM of every 
walker back to the origin. However, the importance-sampled 
propagator involves ratios of overlaps with the trial wave func- 
tion (^rpl^), which would have to be corrected in the random 
walk whenever a shift is made. 

Instead our solution to this diffusive motion is to let the trial 
wave function slide along with the walkers. In other words, 
we rewrite the kinetic energy operator as 

f = f cm +f, (Bl) 

where T cm represents the CM kinetic energy, and T' the inter- 
nal kinetic energy in the CM frame. The total Hamiltonian is 
given by 

H = f cm + f ' + V = T cm + H' . (B2) 

The quantities that we wish to compute are governed by the 
"internal" Hamiltonian H' . Since V involves only relative 
coordinates among the particles, it commutes with T cm ; or 
more generally, 

[f cm , H'] = 0. (B3) 

In this way, the importance-sampled QMC propagation is de- 
termined by H'. The motion of the CM in each walker is 
a separate free diffusion which is governed by T cm . In the 



random-walk process, we are now free to correct for the CM 
motion by shifting the walkers back to the origin whenever 
necessary. For consistency, this correction must be applied 
both in the normal random walk and in the back-propagation 
phase. 

2. Separating the center-of-mass kinetic energy 

The moving trial wave function, however, poses a problem 
for the calculation of the kinetic energy. Now the orbitals are 
free to slide, and the diffusive motion of the orbital's CM is 
no longer suppressed in the LAB frame. When we use the 
usual f-term in the Hamiltonian in Eq. J43l > to compute the ki- 
netic energy, we obtain the total (T), in which T cm = (T cm ) 
and the desired (T") are mixed. This leads to a spurious in- 
crease in the estimate of the kinetic energy and consequently 
the total energy. For example, the uncorrected ground-state 
energy for the system shown in Table UTI would be —1.887(2) 
with (T) = 2.092(3); thus the total energy is overestimated 
by 0.08 due to the contribution from T cm . Since we know the 
nature of the CM motion, it is fairly straightforward to extract 
T cm and explicitly subtract it from the kinetic and total energy 
estimates. Allowing the droplet to freely slide in the calcula- 
tion is equivalent to having a spurious "propagator" e~ ArTcm , 
whose effect on the wave function for the CM is described by 
the diffusion equation 

_«!)=f * (R T ) 

T 

It is a well known property of such a diffusion process that the 
averaged squared distance (R 2 (r)) grows linearly with the 
(imaginary) time r: 

(R 2 (r)) = 6r. 

We can obtain b by recording the quantity (R 2 (t)) for a pe- 
riod of time in the QMC simulation. The constant b is lin- 
early proportional to T cm . More specifically, the center-of- 
mass Hubbard hopping parameter t CTa can be extracted from 
6: 

t cm = b/2 . (B4) 

This gives us the correct kinetic and total energies without the 
spurious center-of-mass motion: 

(f ') = (1 - V) CO i ( B5a > 
(H') = (f) + (V 2B ) . (B5b) 

To conclude, there are two necessary modifications in the 
QMC algorithm in order to treat quantum droplets which are 
not confined: 

1. We let the trial wave function effectively "follow" the 
QMC orbitals, by defining its CM with that of each 
QMC orbital. 
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2. For each orbital, we keep track and accumulate all the 
applied CM shifts in order to estimate (R 2 (r)). This 
gives us the fraction of CM kinetic energy through the 
constant i cm . 



These modifications in the QMC allows us to obtain the cor- 
rect density profile and energies of a translationally-invariant 
Hamiltonian. 
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